System size expansion for systems with an absorbing state 
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The well known van Kampen system size expansion, while of rather general applicability, is 
shown to fail to reproduce some qualitative features of the time evolution for systems with an 
absorbing state, apart from a transient initial time interval. We generalize the van Kampen ansatz 
by introducing a new prescription leading to non-Gaussian fluctuations around the absorbing state. 
The two expansion predictions are explicitly compared for the infinite range voter model with 
speciation as a paradigmatic model with an absorbing state. The new expansion, both for a finite 
size system in the large time limit and at finite time in the large size limit, converges to to the exact 
solution as obtained in a numerical implementation using the Gillespie algorithm. Furthermore, the 
predicted lifetime distribution is shown to have the correct asymptotic behavior. 
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The time evolution of systems consisting of large num- 
ber of discrete entities such as photons, nuclei, proteins 
or organisms is often described by a master equation, 
a differential equation which, in most cases, cannot be 
solved analytically. The van Kampen system-size ex- 
pansion 0, Q is one of the techniques typically used 
to overcome such a limitation, although alternative ap- 
proaches have been proposed |3| ■ This method allows one 
to account for the essential aspects of the problem and 
provides a very useful tool to approximate the temporal 
evolution. However, such an approach is able to charac- 
terize the fluctuations properly when the system has no 
boundaries or evolves far from them Q. For instance, if 
a system is driven towards ultimate extinction, the van 
Kampen approximation is at best appropriate at short 
times. 

When a system with no boundaries initially has a large 
number of particles, one expects that the macroscopic 
evolution is relatively less affected by fluctuations at least 
within a finite temporal scale. This general consideration 
leads to the rule of thumb that deviations from the col- 
lective behavior are of order VN, where N is the size of 
the system. More specifically, the population of the sys- 
tem, n, can be split into two contributions: a macroscopic 
part of order N, Ncj>(t), whose evolution is determinis- 
tic; and a random variable of order v^V, V^V£- This is 
the celebrated van Kampen ansatz, n = N<f>(t) + VNS,, 
which approximates random jumps around the macro- 
scopic part with Gaussian fluctuations and naturally in- 
troduces a small parameter for large N, 1/y/N, that can 
be used as an expansion parameter for the solution of the 



master equation. However, if the system has an absorb- 
ing state it will be driven towards a final absorption, for 
example, eventual extinction. Thus, sooner or later, the 
fluctuations could become comparable with the macro- 
scopic part, despite starting off with a large number of 
individuals. This means that the validity of the van Kam- 
pen approximation may be limited to a short initial time 
interval and fluctuations may no longer be Gaussian. 

As a paradigmatic example of a system with an ab- 
sorbing state, we will consider the infinite range voter 
model with speciation, a simple model that can be han- 
dled analytically. By exploiting standard methods used 
for diffusion processes with absorbing boundaries, we will 
consider an improvement of the classical van Kampen 
technique. However, we will show that despite the mod- 
ification, the new version fails to match numerical simu- 
lations, thus calling for a different approach. To provide 
a general context, we consider the following birth and 
death master equation which is commonly encountered 
in population ecology, 



-P n (t) = (e- 1 -l)[T(n + l\n)P n (t)) 

+ (e+ 1 -l)[T(n-l\n)P n (t)} (1) 

where P n (t) is the probability of observing the system in 
the state n at time t, and the shift operators e^ 1 act on 
a function as e^ 1 /(n) = /(n± 1). T(i\j) is the transition 
probability from state j to state i. In ecology, the state n 
would correspond to a species abundance n. When pop- 
ulations are large, customarily birth and death transition 
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rates turn out to be analytic functions of the density n/N, 
namely T(n±l|n) = T±(n/N), where N denotes the total 
number of individuals or particles. This naturally sug- 
gests a parameter which fluctuations can be compared to. 
Thus, in order to obtain the correct system-size expan- 
sion we introduce the following generalized van Kampen 
ansatz 



n = N<t>{t) + N a £ 



(2) 



with ^ a < 1. Under this assumption n/N = <j>(t) + 
N a ~ 1 t; and when N ^> 1 the transition rates can be 
expanded as power series of y = £iV Q-1 according to 



T±{cj> + y) 



fc=0 



(3) 



where T±\z) is the fc-th derivative of T±(z). Since N is 
large the shift operators have the following representation 
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Substituting Eqs. ©, © and g]) into Eq.((T]) and 
defining the new probability distribution II as II(£, t) oc 
P n {t), one can collect terms proportional to different 
powers of N. In the limit of large N the leading or- 
der provides the usual macroscopic law defined by the 
following deterministic equation 



— = T + (0)-r_(0) 



(5) 



where r = t/N and we assume that T+(0) — T_(<£) is not 
identically zero. Working out the next -to-leading orders, 
one eventually obtains a differential equation which up to 
the second derivative reads: 
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where the time dependence of 4> is given by Eq. ([3]). The 
right-hand side of this equation contains two series which 
are negligible with respect to the first three terms when 
N is large which are proportional to N°, N 1 ^ 2 " and N~ a 



respectively. If we assume that both T ( y{<f)-T ( +> (</>) and 
T—(<j)) + ?+(</>) are different from zero in order to avoid 
the trivial result of vanishing fluctuations, one has to set 



a = 1/2 in Eq.© Accordingly, in the limit N 
recover the standard van Kampen equation 



oo we 
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T_(4>)+T + (ct>) (7) 



if 



which is a linear Fokkcr-Planck equation whose solution 
is a non-stationary Gaussian distribution. 

However, for systems with absorbing boundaries at 
n = and large temporal scales, the term proportional 
to iV 1-2 ", T-{4>) + T+(</>), approaches zero while the one 
proportional to N~ a , T { _}\(j>) + (</>), does not. This 
is because when 0—^0, ± T + ((j>) are proportional 

to (f> (in the case of a simple absorbing state). In this 
case the van Kampen prescription is no longer valid and 
we need to set a = in the limit of large N. For these 
systems fluctuations are progressively more important in 
the long run, because N(f>(r) -c £. Thus, the differential 
equation governing the fluctuations is well approximated 
by the following Fokker-Planck equation 



d_ 



n 



T^>(0) -Tj ij (0) 



(i) 



T (i) (0) + T (i) (0) ]^^ n ) (8) 



This equation is different from Eq. ([7J owing to the linear 
diffusion term which results in the fluctuations being no 
longer Gaussian distributed. Furthermore, if we do not 
provide Eq. ([7]) with an absorbing boundary condition 
for the solution, fluctuations could lead to negative val- 
ues for n. In contrast, Eq. @ has a natural boundary at 
£ = which prevents fluctuations, and thus n, from be- 
coming negative. Similar considerations also hold when 
transition rates have more general alg ebraic behaviors in 
the vicinity of the absorbing state [15| . Interestingly, Eq. 
(HJ) can be exactly solved @, its solution being 



n(£,r|£ o ,0) = 



1 



D 1 



where fx = T_ (0) — T+'(0) is supposed to be positive, 

D = [tP(0) + T} 1) (0)]/2, I x (z) is the modified Bessel 
function of the first kind and £o is the value of £ when 
t = 0. It is worth noting that although the solution is 
absorbing, one gets Iim^o fl(£, r|£ , 0) ^ 0. 
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The theory we have developed so far can be applied 
straightforwardly to many different absorbing systems. 
In particular, we now focus on the infinite range voter 
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FIG. 1: Probability distribution for the fluctuations £ at dif- 
ferent times and with different parameters specified in the 
top right of the panels. The noisy lines are the results of nu- 
merical solutions obtained by averaging over 10 6 stochastic 
realizations, the solid-dotted lines represent the Gaussian so- 
lution modified for absorbing boundaries given by (|16[) . Fluc- 
tuations are approximately Gaussian distributed only for rel- 
atively short times, while for large r the Gaussian assumption 
breaks down. In all cases v = O.Of and no = 300. 



model with speciation, a particular case of the more gen- 
eral voter model which is of interest in o pini on formation 
problems |^p, but also in biological [9l.ll(| and ecologi- 
cal contexts (ill . \v\ . 

The modified version of the voter model we investigate 
is characterized by a parameter, the speciation rate v, 
which averts the collapse of the whole system into a triv- 
ial monodominant state characterized by <fi = 1. Specifi- 
cally, let us consider a system composed of TV elements, 
all of them mutually interacting and belonging to possi- 
bly different species. If we now focus on a specific species, 
we can re-map all elements with two labels: the label X\ 
for the elements of the selected species, the label Xq for 
the rest. Finally, at each time step we randomly choose 
and update a pair of elements according to the following 
interaction rules: 



X\ + Xq > Xq + Xq 

Xq + Xi > X\ + X\ 

X\ + X\ > X\ + Xq. 



(10) 

(11) 

(12) 



An individual of the species of the first term on the lhs is 
envisaged to be replaced by an individual of the second 
term on the lhs except for speciation which occurs with 
a probability v as in the third rule. The factors above 
the arrow denote the probability of the event indicated 
in the equation. 

Let us denote by n the number of X\ individuals, so 
that TV — n is the total number of elements of type Xo. 
According to (|10p - (fT2"l) the only transitions allowed are 
those from n to n =b 1, and the corresponding transition 



probabilities read 
T(n — l|n) 
T(n +l|n) 



n N — n n 
{1 - V) NN—1 +V N 



TV — n n 
TV TV- 1' 



(13) 
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where the initial states are on the right and final states 
on the left. Since T(±1|0) = 0, once the population of 
Xi dies out, the selected species cannot be re-introduced 
into the system. Thus, this model has a continual turn 
over of species: new species appear at rate v, but eventu- 
ally they go extinct. This implies that n = is an absorb- 
ing state. On the contrary, when the population of X\ 
reaches the maximum value TV, transitions to TV + 1 are 
not allowed since T(TV + 1|TV) = 0, while T(TV-1|TV) = v. 
As a consequence, n = N is a reflecting boundary. In the 
following wc will focus on the time evolution of the sys- 
tem when n is kept finite as TV becomes larger and larger 

M- 

If we apply the generalized expansion described in the 
previous section, we find that the macroscopic law ac- 
cording to Eq. ([5]) is 4> = —vcf), thus (T>(t) = 4>oe~ l/T with 
t = t/(N — 1). The van Kampcn equation corresponding 
to Eq. reads 
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where /(t) = 1/2 [(2 
Its solution is 
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with v (t) = (2 - is)Me ur - l)/(2f) - (1 - vWqT. 

In order to account for the absorbing boundary, we 
added a time dependent constraint on £, so that n in 
Eq. @ varies between and TV. To guarantee this 
latter condition, wc imposed £ m i„ ^ £ $5 ( rara , where 
£mm = -VN<t> e-" T and Uax = VN(1- 4> e- VT ). In 
correspondence to n = 0, £ = £ m ;„ is an absorbing 
boundary, while £ = £maxi which corresponds to n = TV, 
is a reflecting boundary. The final solution accounting 
for the absorbing boundary reads 



IW£,t) 



\/47r 7/(t) 



exp 
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This solution has a delta peak at £ = (or n 
t — y and vanishes at £ = — \/TV0o e ~ VT ■ 
We now turn to Eq. ((5J which now reads 
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FIG. 2: Comparison between the numerical simulations of 
the rules in Eqs. (|10[1 - H12[) and the theoretical predictions 
of the probability distribution in Eq. @ at different times 
and systems sizes, according to the parameters specified in 
the top right of the panels. The noisy lines are the numerical 
profiles obtained by averaging over f 6 stochastic simulations. 
While initially fluctuations are approximately Gaussian, in 
the sequel they are non-Gaussian and peak at zero. In all 
cases v — 0.01 and no = 300. 



where we have used that T_ (0) — T]_ l; (0) = v and 



T} lj (0) = 2-v (see Eqs. (O and (fli])) 



T«(0) 

In order to test the validity of the methods, we per- 
formed extensive numerical simulations of Eqs. (jTT))) - (fT2")) 
through the Gillespie algorithm [l3| , which allows one to 
produce time series which exactly recover the solution 
of the master equation ^ with the rates in Eqs. (fT3"|) 
and (Q3J • Fig. ([1]) shows typical results of the stochas- 
tic simulations and their comparison with the absorbing 
van Kampen solution in Eq. (|16[) . For short times the 
first three profiles overlap well, but as time increases, 
the probability distribution does not match the numer- 
ical simulation, as shown in the last three panels. Fur- 
thermore, the agreement does not improve on increasing 
the size of the system. In contrast, Fig. ^ shows the 
solution in Eq. © with fj, = v and D = (2 — v)/2. In- 



creasing TV, while keeping r fixed, improves the matching 
as shown in the first three panels. As expected Eq. © 
converges to the numerical profiles as r increases. 



As illustrated in both figures, in the presence of an 
absorbing state, the system is characterized by at least 
two temporal scales, n and T2, which make fluctua- 
tions evolve according to Eq. ([7]) for r < T\ and Eq. 
(O for r > T2- It is possible to estimate roughly 
the two scales by observing that one should expect 
the generalized van Kampen ansatz to work until the 
fluctuations are of the same magnitude as the macro- 
scopic part, namely N(f)(r) ~ N a a(r) where <j(t) is 
the variance of £. When a = 1/2, this condition 
translates into N(j)Qe~ VT ~ y / 2A?7(r) which gives n ~ 
l/(3f) ln(not'/(2— u)) fori/r ^> 1. For the expansion with 
a = 0, we have N^ Q e~ UT ^ e _J/T - v /(2 - v)n /v{e VT - 1) 
which gives — l/v\n(i/no/(2 — v) + 1). Note that 
Ti < T2- The above mentioned condition of validity of 
the classic van Kampen expansion is confirmed by nu- 
merical simulations (data not shown). 

Finally, for systems with absorbing boundaries, it is in- 
teresting to calculate an analytical expression for the sur- 
vival probability Ps(t) [Oj]- I n our case, we get the exact 
expression P s (t) = 1 -'[(1 - e~ VT )/(l - (1 - v)e- VT )] n ° 
which, in the scaling limit v — > with vt fixed, simpli- 
fies to f{vr)/t with f(z) = zn /(e z — 1). Amazingly, 
the same result is also obtainable using Eq. ©, with 
P s (r)=/den(£,r|6,,0). 

Summarizing, in the presence of systems with absorb- 
ing states, one has to generalize the standard van Kam- 
pen ansatz in order to monitor the temporal evolution at 
large times. As time elapses, fluctuations become more 
and more important and are no longer Gaussian. How- 
ever, they still can be analytically treated and lead to the 
general solution given by Eq. (|9]). 
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